Goal: can machine learning methods help us to associate metabolites with leaf length?

Previously (script 11b2) I filtered out unnamed metabolites. Here I keep them all.

Also I will PC separately for root and leaf.

library(glmnet)
Loading required package: Matrix
Loaded glmnet 4.0-2
library(relaimpo)
Loading required package: MASS
Loading required package: boot
Loading required package: survey
Loading required package: grid
Loading required package: survival

Attaching package: ‘survival’

The following object is masked from ‘package:boot’:

    aml


Attaching package: ‘survey’

The following object is masked from ‘package:graphics’:

    dotchart

Loading required package: mitools
This is the global version of package relaimpo.

If you are a non-US user, a version with the interesting additional metric pmvd is available

from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.3     ✓ purrr   0.3.4
✓ tibble  3.0.4     ✓ dplyr   1.0.2
✓ tidyr   1.1.2     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x tidyr::expand() masks Matrix::expand()
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x tidyr::pack()   masks Matrix::pack()
x dplyr::select() masks MASS::select()
x tidyr::unpack() masks Matrix::unpack()

get leaflength data

leaflength <- read_csv("../../plant/output/leaf_lengths_metabolite.csv") %>%
  mutate(pot=str_pad(pot, width=3, pad="0"),
         sampleID=str_c("wyo", genotype, pot, sep="_")) %>%
  select(sampleID, leaf_avg_std) 

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  pot = col_double(),
  soil = col_character(),
  genotype = col_character(),
  trt = col_character(),
  leaf_avg = col_double(),
  leaf_avg_std = col_double()
)
leaflength %>% arrange(sampleID)

get and wrangle metabolite data

met_raw <-read_csv("../input/metabolites_set1.csv")

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_double(),
  tissue = col_character(),
  soil = col_character(),
  genotype = col_character(),
  autoclave = col_character(),
  time_point = col_character(),
  concatenate = col_character()
)
ℹ Use `spec()` for the full column specifications.
met <- met_raw %>% 
  mutate(pot=str_pad(pot, width = 3, pad = "0")) %>%
  mutate(sampleID=str_c("wyo", genotype, pot, sep="_")) %>%
  select(sampleID, genotype, tissue, sample_mass = `sample_mass mg`, !submission_number:concatenate) %>%
  pivot_longer(!sampleID:sample_mass, names_to = "metabolite", values_to = "met_amount") %>%
  
  #adjust by sample mass
  mutate(met_per_mg=met_amount/sample_mass) %>%
  
  #scale and center
  group_by(metabolite, genotype, tissue) %>%
  mutate(met_per_mg=scale(met_per_mg),
         met_amt=scale(met_amount)
  ) %>% 
  pivot_wider(id_cols = sampleID, 
              names_from = c(tissue, metabolite), 
              values_from = starts_with("met_"),
              names_sep = "_")

met 

split this into two data frames, one normalized by tissue amount and one not.

met_per_mg <- met %>% select(sampleID,  starts_with("met_per_mg")) %>%
  as.data.frame() %>% column_to_rownames("sampleID")
met_amt <- met %>% select(sampleID,  starts_with("met_amt")) %>%
  as.data.frame() %>% column_to_rownames("sampleID")

get leaf data order to match

leaflength <- leaflength[match(met$sampleID, leaflength$sampleID),]
leaflength

Calc PCAs:

normalized

leaf

met_per_mg.leaf_PCA <- met_per_mg %>% 
  select(matches("_leaf_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.leaf_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_per_mg.leaf_PCA$sdev^2, PC=str_c("PC", 
                                                      str_pad(1:length(met_per_mg.leaf_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, normalized leaf metabolites")

root

met_per_mg.root_PCA <- met_per_mg %>% 
  select(matches("_root_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.root_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_per_mg.root_PCA$sdev^2, PC=str_c("PC", 
                                                      str_pad(1:length(met_per_mg.root_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, normalized root metabolites")

raw

leaf

met_amt.leaf_PCA <- met_amt %>%
  select(matches("_leaf_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.leaf_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_amt.leaf_PCA$sdev^2, PC=str_c("PC", 
                                                   str_pad(1:length(met_amt.leaf_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, raw leaf metabolites")

root

met_amt.root_PCA <- met_amt %>%
  select(matches("_root_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.root_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_amt.root_PCA$sdev^2, PC=str_c("PC", 
                                                   str_pad(1:length(met_amt.root_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, raw root metabolites")

now try these in a penalized regression

normalized

are the PCs normalized?

colMeans(met_amt.leaf_PCA$x) %>% round(3) #yes centered
 PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 
   0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
PC27 PC28 PC29 PC30 PC31 PC32 PC33 PC34 PC35 PC36 
   0    0    0    0    0    0    0    0    0    0 
apply(met_amt.leaf_PCA$x, 2, sd) %>% round(2) #not scaled
  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13  PC14  PC15  PC16  PC17  PC18  PC19  PC20  PC21  PC22 
11.94  8.37  6.74  5.84  5.46  5.32  5.02  4.55  4.38  4.11  3.87  3.80  3.73  3.56  3.44  3.32  3.27  3.17  3.16  3.08  2.94  2.90 
 PC23  PC24  PC25  PC26  PC27  PC28  PC29  PC30  PC31  PC32  PC33  PC34  PC35  PC36 
 2.86  2.78  2.73  2.61  2.56  2.53  2.40  2.39  2.33  2.26  2.22  2.14  0.00  0.00 

combine the leaf and root, and then scale them:

met_per_mg.leaf_PCs <- met_per_mg.leaf_PCA$x
colnames(met_per_mg.leaf_PCs) <- str_c("leaf_", colnames(met_per_mg.leaf_PCs))

met_per_mg.root_PCs <- met_per_mg.root_PCA$x
colnames(met_per_mg.root_PCs) <- str_c("root_", colnames(met_per_mg.root_PCs))

met_per_mg.PCs <- cbind(met_per_mg.leaf_PCs, met_per_mg.root_PCs) %>%
  scale()

met_amt.leaf_PCs <- met_amt.leaf_PCA$x
colnames(met_amt.leaf_PCs) <- str_c("leaf_", colnames(met_amt.leaf_PCs))

met_amt.root_PCs <- met_amt.root_PCA$x
colnames(met_amt.root_PCs) <- str_c("root_", colnames(met_amt.root_PCs))

met_amt.PCs <- cbind(met_amt.leaf_PCs, met_amt.root_PCs) %>%
  scale()

also combine the rotations

met_per_mg.leaf_rotation <- met_per_mg.leaf_PCA$rotation %>%
  as.data.frame() %>% 
  rename_with(~ str_c("leaf_", .x)) %>%
  rownames_to_column("metabolite")

met_per_mg.root_rotation <- met_per_mg.root_PCA$rotation %>%
  as.data.frame() %>% 
  rename_with(~ str_c("root_", .x)) %>%
  rownames_to_column("metabolite")

met_per_mg.PC_rotation <- full_join(met_per_mg.leaf_rotation, met_per_mg.root_rotation, by="metabolite")

met_amt.leaf_rotation <- met_amt.leaf_PCA$rotation %>% 
  as.data.frame() %>% 
  rename_with(~ str_c("leaf_", .x)) %>%
  rownames_to_column("metabolite")

met_amt.root_rotation <- met_amt.root_PCA$rotation %>%
  as.data.frame() %>% 
  rename_with(~ str_c("root_", .x)) %>%
  rownames_to_column("metabolite")

met_amt.PC_rotation <- full_join(met_amt.leaf_rotation, met_amt.root_rotation, by="metabolite")
met_per_mg_fit1LOO <- cv.glmnet(x=met_per_mg.PCs, y=leaflength$leaf_avg_std, nfolds = nrow(met_per_mg.PCs), alpha=1 )
Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold
plot(met_per_mg_fit1LOO)

bestlam=met_per_mg_fit1LOO$lambda.1se

NEXT STEP: Do a K-fold CV, repeat many times and average. Might as well do alpha while we are at it. If we are doing alpha, then we need to manually create our own folds list for each run

normalized

multi CV

Fit 101 CVs for each of 11 alphas

set.seed(1245)

folds <- tibble(run=1:101) %>% 
  mutate(folds=map(run, ~ sample(rep(1:6,6))))

system.time (met_per_mg_multiCV <- expand_grid(run=1:100, alpha=round(seq(0,1,.1),1)) %>%
               left_join(folds, by="run") %>%
               mutate(fit=map2(folds, alpha, ~ cv.glmnet(x=met_per_mg.PCs, y=leaflength$leaf_avg_std, foldid = .x, alpha=.y
                                                         )))
             #, lambda=exp(seq(-5,0,length.out = 50)) )))
) #100 seconds
   user  system elapsed 
 93.102   4.877 118.425 
head(met_per_mg_multiCV)

for each fit, pull out the mean cv error, lambda, min lambda, and 1se lambda

met_per_mg_multiCV <- met_per_mg_multiCV %>%
  mutate(cvm=map(fit, magrittr::extract("cvm")),
         lambda=map(fit, magrittr::extract("lambda")),
         lambda.min=map_dbl(fit, magrittr::extract("lambda.min" )),
         lambda.1se=map_dbl(fit, magrittr::extract("lambda.1se")),
         nzero=map(fit, magrittr::extract("nzero"))
  )

head(met_per_mg_multiCV)

now calculate the mean and sem of cvm and min,1se labmdas. These need to be done separately because of the way the grouping works

met_per_mg_summary_cvm <- met_per_mg_multiCV %>% dplyr::select(-fit, -folds) %>% 
  unnest(c(cvm, lambda)) %>%
  group_by(alpha, lambda) %>%
  summarize(meancvm=mean(cvm), sem=sd(cvm)/sqrt(n()), high=meancvm+sem, low=meancvm-sem)
`summarise()` regrouping output by 'alpha' (override with `.groups` argument)
met_per_mg_summary_cvm
met_per_mg_summary_lambda <- met_per_mg_multiCV %>% dplyr::select(-fit, -folds, -cvm) %>% 
  group_by(alpha) %>%
  summarize(
    lambda.min.sd=sd(lambda.min), 
    lambda.min.mean=mean(lambda.min),
    #lambda.min.med=median(lambda.min), 
    lambda.min.high=lambda.min.mean+lambda.min.sd,
    #lambda.min.low=lambda.min.mean-lambda.min.sem,
    #lambda.1se.sem=sd(lambda.1se)/sqrt(n()), 
    lambda.1se.mean=mean(lambda.1se),
    #lambda.1se.med=median(lambda.1se), 
    #lambda.1se.high=lambda.1se+lambda.1se.sem,
    #lambda.1se.low=lambda.1se-lambda.1se.sem,
    nzero=nzero[1],
    lambda=lambda[1]
  )
`summarise()` ungrouping output (override with `.groups` argument)
met_per_mg_summary_lambda

plot it

met_per_mg_summary_cvm %>%
  #filter(alpha!=0) %>% # worse than everything else and throwing the plots off
  ggplot(aes(x=log(lambda), y= meancvm,  ymin=low, ymax=high)) +
  geom_ribbon(alpha=.25) +
  geom_line(aes(color=as.character(alpha))) +
  facet_wrap(~ as.character(alpha)) +
   coord_cartesian(xlim=(c(-5,0))) +
  geom_vline(aes(xintercept=log(lambda.min.mean)), alpha=.5, data=met_per_mg_summary_lambda) +
  geom_vline(aes(xintercept=log(lambda.min.high)), alpha=.5, data=met_per_mg_summary_lambda, color="blue") 

So overall these look more reasonable than the LOO plot.

Make a plot of MSE at minimum lambda for each alpha

met_per_mg_summary_cvm %>% 
  group_by(alpha) %>%
  filter(rank(meancvm, ties.method = "first")==1) %>%
  ggplot(aes(x=alpha,y=meancvm,ymin=low,ymax=high)) +
  geom_ribbon(color=NA, fill="gray80") +
  geom_line() +
  geom_point()

not a particular large difference there, aside from 0.1 and even then, not too much better.

Plot the number of nzero coefficients

met_per_mg_summary_lambda %>%
  unnest(c(lambda, nzero)) %>%
  group_by(alpha) %>%
  filter(abs(lambda.min.mean-lambda)==min(abs(lambda.min.mean-lambda))  ) %>%
  ungroup() %>%

ggplot(aes(x=as.character(alpha), y=nzero)) +
  geom_point() +
  ggtitle("Number of non-zero coefficents at minimum lambda") +
  ylim(0,36)

OK let’s do repeated test train starting from these CV lambdas

multi_tt <- function(lambda, alpha, n=10000, sample_size=36, train_size=30, x, y=leaflength$leaf_avg_std) {
  print(lambda)
  print(alpha)
tt <-
  tibble(run=1:n) %>%
  mutate(train=map(run, ~ sample(1:sample_size, train_size))) %>%
  mutate(fit=map(train, ~ glmnet(x=x[.,], y=y[.], lambda = lambda, alpha = alpha ))) %>%
  
  mutate(pred=map2(fit, train, ~ predict(.x, newx = x[-.y,]))) %>%
  mutate(cor=map2_dbl(pred, train, ~ cor(.x, y[-.y])  )) %>%
  mutate(MSE=map2_dbl(pred, train, ~ mean((y[-.y] - .x)^2))) %>%
  summarize(
    num_na=sum(is.na(cor)), 
    num_lt_0=sum(cor<=0, na.rm=TRUE),
    avg_cor=mean(cor, na.rm=TRUE),
    avg_MSE=mean(MSE))
tt
}

per_mg_fit_test_train <- met_per_mg_summary_lambda %>% 
  select(alpha, lambda.min.mean)

per_mg_fit_test_train <- met_per_mg_multiCV %>%
  filter(run==1) %>%
  select(alpha, fit) %>%
  right_join(per_mg_fit_test_train)
Joining, by = "alpha"
per_mg_fit_test_train <- per_mg_fit_test_train %>%
  mutate(pred_full=map2(fit, lambda.min.mean, ~ predict(.x, s=.y, newx=met_per_mg.PCs)),
         full_R=map_dbl(pred_full, ~ cor(.x, leaflength$leaf_avg_std)),
         full_MSE=map_dbl(pred_full, ~ mean((leaflength$leaf_avg_std-.x)^2))) %>%
  
  mutate(tt=map2(lambda.min.mean, alpha, ~ multi_tt(lambda=.x, alpha=.y, x=met_per_mg.PCs)))
[1] 13.36988
[1] 0
[1] 0.7666849
[1] 0.1
[1] 0.8732274
[1] 0.2
[1] 0.6950232
[1] 0.3
[1] 0.6122284
[1] 0.4
[1] 0.5404081
[1] 0.5
[1] 0.4719584
[1] 0.6
[1] 0.4194147
[1] 0.7
[1] 0.3761253
[1] 0.8
[1] 0.3410749
[1] 0.9
[1] 0.3097519
[1] 1
(per_mg_fit_test_train <- per_mg_fit_test_train %>% unnest(tt))
per_mg_fit_test_train %>%
  ggplot(aes(x=alpha)) +
  geom_line(aes(y=avg_cor), color="red") +
  geom_point(aes(y=avg_cor), color="red") +
  geom_line(aes(y=avg_MSE), color="blue") +
  geom_point(aes(y=avg_MSE), color="blue")

alpha of 0.8 to 1.0 are very similar and are the best here.

look at fit:

alpha_per_mg <- .8

best_per_mg <- per_mg_fit_test_train %>% filter(alpha == alpha_per_mg) 
best_per_mg_fit <- best_per_mg$fit[[1]]
best_per_mg_lambda <- best_per_mg$lambda.min.mean

per_mg_coef.tb <- coef(best_per_mg_fit, s=best_per_mg_lambda) %>% 
  as.matrix() %>% as.data.frame() %>% 
  rownames_to_column(var="PC") %>%
  rename(beta=`1`)
  
per_mg_coef.tb %>% filter(beta!=0) %>% arrange(beta)
NA

pred and obs

plot(leaflength$leaf_avg_std, best_per_mg$pred_full[[1]])

cor.test(leaflength$leaf_avg_std, best_per_mg$pred_full[[1]]) #.57

    Pearson's product-moment correlation

data:  leaflength$leaf_avg_std and best_per_mg$pred_full[[1]]
t = 4.1276, df = 34, p-value = 0.0002243
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3076242 0.7617163
sample estimates:
      cor 
0.5777675 
best_per_mg$full_MSE
[1] 0.7717674

Percent variance explained

per_mg_vars <- per_mg_coef.tb %>% 
  filter(beta !=0, PC!="(Intercept)") %>%
  pull(PC) %>% c("leaf_avg_std", .)

per_mg_relimp <- leaflength %>% select(leaf_avg_std) %>% cbind(met_per_mg.PCs) %>% as.data.frame() %>% dplyr::select(all_of(per_mg_vars)) %>%
  calc.relimp() 

per_mg_coef.tb <- per_mg_relimp@lmg %>% as.matrix() %>% as.data.frame() %>%
  rownames_to_column("PC") %>%
  rename(PropVar_met_per_mg=V1) %>%
  full_join(per_mg_coef.tb) %>%
  arrange(desc(PropVar_met_per_mg))
Joining, by = "PC"
per_mg_coef.tb
NA

Checkout the rotations.

met_per_mg_rotation_out <- met_per_mg.PC_rotation %>% 
  pivot_longer(-metabolite, names_to="PC", values_to="loading") %>%
  filter(PC %in% filter(per_mg_coef.tb, beta!=0)$PC ) %>%
  group_by(PC) %>%
    filter(!str_detect(metabolite,".*(leaf|root)_[0-9]*$")) %>%
  filter(abs(loading) >= 0.05) %>%
  left_join(per_mg_coef.tb, by="PC") %>%
  arrange(desc(abs(beta)), desc(abs(loading))) %>%
  mutate(organ=ifelse(str_detect(metabolite, "_leaf_"), "leaf", "root"),
         transformation="normalized",
         metabolite=str_remove(metabolite, "met_per_mg_(root|leaf)_"),
         metabolite_effect_on_leaf=ifelse(beta*loading>0, "increase", "decrease"))
met_per_mg_rotation_out %>%  write_csv("../output/Leaf_associated_metabolites_normalized.csv")

met_per_mg_rotation_out

non-normazlized

multi CV

Fit 101 CVs for each of 11 alphas

set.seed(1245)

folds <- tibble(run=1:101) %>% 
  mutate(folds=map(run, ~ sample(rep(1:6,6))))

system.time (met_amt_multiCV <- expand_grid(run=1:100, alpha=round(seq(0,1,.1),1)) %>%
               left_join(folds, by="run") %>%
               mutate(fit=map2(folds, alpha, ~ cv.glmnet(x=met_amt.PCs, y=leaflength$leaf_avg_std, foldid = .x, alpha=.y
                                                         )))
             #, lambda=exp(seq(-5,0,length.out = 50)) )))
) #100 seconds
   user  system elapsed 
 70.687   2.897  87.024 
head(met_amt_multiCV)

for each fit, pull out the mean cv error, lambda, min lambda, and 1se lambda

met_amt_multiCV <- met_amt_multiCV %>%
  mutate(cvm=map(fit, magrittr::extract("cvm")),
         lambda=map(fit, magrittr::extract("lambda")),
         lambda.min=map_dbl(fit, magrittr::extract("lambda.min" )),
         lambda.1se=map_dbl(fit, magrittr::extract("lambda.1se")),
         nzero=map(fit, magrittr::extract("nzero"))
  )

head(met_amt_multiCV)

now calculate the mean and sem of cvm and min,1se labmdas. These need to be done separately because of the way the grouping works

met_amt_summary_cvm <- met_amt_multiCV %>% dplyr::select(-fit, -folds) %>% 
  unnest(c(cvm, lambda)) %>%
  group_by(alpha, lambda) %>%
  summarize(meancvm=mean(cvm), sem=sd(cvm)/sqrt(n()), high=meancvm+sem, low=meancvm-sem)
`summarise()` regrouping output by 'alpha' (override with `.groups` argument)
met_amt_summary_cvm
met_amt_summary_lambda <- met_amt_multiCV %>% dplyr::select(-fit, -folds, -cvm) %>% 
  group_by(alpha) %>%
  summarize(
    lambda.min.sd=sd(lambda.min), 
    lambda.min.mean=mean(lambda.min),
    #lambda.min.med=median(lambda.min), 
    lambda.min.high=lambda.min.mean+lambda.min.sd,
    #lambda.min.low=lambda.min.mean-lambda.min.sem,
    #lambda.1se.sem=sd(lambda.1se)/sqrt(n()), 
    lambda.1se.mean=mean(lambda.1se),
    #lambda.1se.med=median(lambda.1se), 
    #lambda.1se.high=lambda.1se+lambda.1se.sem,
    #lambda.1se.low=lambda.1se-lambda.1se.sem,
    nzero=nzero[1],
    lambda=lambda[1]
  )
`summarise()` ungrouping output (override with `.groups` argument)
met_amt_summary_lambda

plot it

met_amt_summary_cvm %>%
  #filter(alpha!=0) %>% # worse than everything else and throwing the plots off
  ggplot(aes(x=log(lambda), y= meancvm,  ymin=low, ymax=high)) +
  geom_ribbon(alpha=.25) +
  geom_line(aes(color=as.character(alpha))) +
  facet_wrap(~ as.character(alpha)) +
   coord_cartesian(xlim=(c(-5,0))) +
  geom_vline(aes(xintercept=log(lambda.min.mean)), alpha=.5, data=met_amt_summary_lambda) +
  geom_vline(aes(xintercept=log(lambda.min.high)), alpha=.5, data=met_amt_summary_lambda, color="blue") 

Make a plot of MSE at minimum lambda for each alpha

met_amt_summary_cvm %>% 
  group_by(alpha) %>%
  filter(rank(meancvm, ties.method = "first")==1) %>%
  ggplot(aes(x=alpha,y=meancvm,ymin=low,ymax=high)) +
  geom_ribbon(color=NA, fill="gray80") +
  geom_line() +
  geom_point()

not a particular large difference here after 0.2

Plot the number of nzero coefficients

met_amt_summary_lambda %>%
  unnest(c(lambda, nzero)) %>%
  group_by(alpha) %>%
  filter(abs(lambda.min.mean-lambda)==min(abs(lambda.min.mean-lambda))  ) %>%
  ungroup() %>%

ggplot(aes(x=as.character(alpha), y=nzero)) +
  geom_point() +
  ggtitle("Number of non-zero coefficents at minimum lambda") +
  ylim(0,36)

OK let’s do repeated test train starting from these CV lambdas

multi_tt <- function(lambda, alpha, n=10000, sample_size=36, train_size=30, x, y=leaflength$leaf_avg_std) {
  print(lambda)
  print(alpha)
tt <-
  tibble(run=1:n) %>%
  mutate(train=map(run, ~ sample(1:sample_size, train_size))) %>%
  mutate(fit=map(train, ~ glmnet(x=x[.,], y=y[.], lambda = lambda, alpha = alpha ))) %>%
  
  mutate(pred=map2(fit, train, ~ predict(.x, newx = x[-.y,]))) %>%
  mutate(cor=map2_dbl(pred, train, ~ cor(.x, y[-.y])  )) %>%
  mutate(MSE=map2_dbl(pred, train, ~ mean((y[-.y] - .x)^2))) %>%
  summarize(
    num_na=sum(is.na(cor)), 
    num_lt_0=sum(cor<=0, na.rm=TRUE),
    avg_cor=mean(cor, na.rm=TRUE),
    avg_MSE=mean(MSE))
tt
}

amt_fit_test_train <- met_amt_summary_lambda %>% 
  select(alpha, lambda.min.mean)

amt_fit_test_train <- met_amt_multiCV %>%
  filter(run==1) %>%
  select(alpha, fit) %>%
  right_join(amt_fit_test_train)
Joining, by = "alpha"
amt_fit_test_train <- amt_fit_test_train %>%
  mutate(pred_full=map2(fit, lambda.min.mean, ~ predict(.x, s=.y, newx=met_amt.PCs)),
         full_R=map_dbl(pred_full, ~ cor(.x, leaflength$leaf_avg_std)),
         full_MSE=map_dbl(pred_full, ~ mean((leaflength$leaf_avg_std-.x)^2))) %>%
  
  mutate(tt=map2(lambda.min.mean, alpha, ~ multi_tt(lambda=.x, alpha=.y, x=met_amt.PCs)))
[1] 121.4399
[1] 0
[1] 0.9724185
[1] 0.1
[1] 0.6990522
[1] 0.2
[1] 0.586571
[1] 0.3
[1] 0.5167138
[1] 0.4
[1] 0.444683
[1] 0.5
[1] 0.3928743
[1] 0.6
[1] 0.3464051
[1] 0.7
[1] 0.3173148
[1] 0.8
[1] 0.2852059
[1] 0.9
[1] 0.2623291
[1] 1
(amt_fit_test_train <- amt_fit_test_train %>% unnest(tt))
amt_fit_test_train %>%
  ggplot(aes(x=alpha)) +
  geom_line(aes(y=avg_cor), color="red") +
  geom_point(aes(y=avg_cor), color="red") +
  geom_line(aes(y=avg_MSE), color="blue") +
  geom_point(aes(y=avg_MSE), color="blue")

alpha of 0.8 to 1.0 are very similar and are the best here.

look at fit:

alpha_amt <- .8

best_amt <- amt_fit_test_train %>% filter(alpha == alpha_amt) 
best_amt_fit <- best_amt$fit[[1]]
best_amt_lambda <- best_amt$lambda.min.mean

amt_coef.tb <- coef(best_amt_fit, s=best_amt_lambda) %>% 
  as.matrix() %>% as.data.frame() %>% 
  rownames_to_column(var="PC") %>%
  rename(beta=`1`)
  
amt_coef.tb %>% filter(beta!=0) %>% arrange(beta)
NA

pred and obs

plot(leaflength$leaf_avg_std, best_amt$pred_full[[1]])

cor.test(leaflength$leaf_avg_std, best_amt$pred_full[[1]]) #.736

    Pearson's product-moment correlation

data:  leaflength$leaf_avg_std and best_amt$pred_full[[1]]
t = 6.3375, df = 34, p-value = 3.15e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5372665 0.8571965
sample estimates:
      cor 
0.7359065 
best_amt$full_MSE
[1] 0.6064862

Percent variance explained

amt_vars <- amt_coef.tb %>% 
  filter(beta !=0, PC!="(Intercept)") %>%
  pull(PC) %>% c("leaf_avg_std", .)

amt_relimp <- leaflength %>% select(leaf_avg_std) %>% cbind(met_amt.PCs) %>% as.data.frame() %>% dplyr::select(all_of(amt_vars)) %>%
  calc.relimp() 

amt_coef.tb <- amt_relimp@lmg %>% as.matrix() %>% as.data.frame() %>%
  rownames_to_column("PC") %>%
  rename(PropVar_met_amt=V1) %>%
  full_join(amt_coef.tb) %>%
  arrange(desc(PropVar_met_amt))
Joining, by = "PC"
amt_coef.tb
NA

Checkout the rotations.

met_amt_rotation_out <- met_amt.PC_rotation %>% 
  pivot_longer(-metabolite, names_to="PC", values_to="loading") %>%
  filter(PC %in% filter(amt_coef.tb, beta!=0)$PC ) %>%
  group_by(PC) %>%
    filter(!str_detect(metabolite,".*(leaf|root)_[0-9]*$")) %>%
  filter(abs(loading) >= 0.05) %>%
  left_join(amt_coef.tb, by="PC") %>%
  arrange(desc(abs(beta)), desc(abs(loading))) %>%
  mutate(organ=ifelse(str_detect(metabolite, "_leaf_"), "leaf", "root"),
         transformation="raw",
         metabolite=str_remove(metabolite, "met_amt_(root|leaf)_"),
         metabolite_effect_on_leaf=ifelse(beta*loading>0, "increase", "decrease"))
met_amt_rotation_out %>%  write_csv("../output/Leaf_associated_metabolites_raw.csv")

met_amt_rotation_out
LS0tCnRpdGxlOiAiTWFjaGluZSBsZWFybmluZyBmb3IgbWV0YWJvbGl0ZXMiCmF1dGhvcjogIkp1bGluIE1hbG9vZiIKZGF0ZTogIjEyLzQvMjAyMCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKYGBgCgpHb2FsOiBjYW4gbWFjaGluZSBsZWFybmluZyBtZXRob2RzIGhlbHAgdXMgdG8gYXNzb2NpYXRlIG1ldGFib2xpdGVzIHdpdGggbGVhZiBsZW5ndGg/CgpQcmV2aW91c2x5IChzY3JpcHQgMTFiMikgSSBmaWx0ZXJlZCBvdXQgdW5uYW1lZCBtZXRhYm9saXRlcy4gIEhlcmUgSSBrZWVwIHRoZW0gYWxsLgoKQWxzbyBJIHdpbGwgUEMgc2VwYXJhdGVseSBmb3Igcm9vdCBhbmQgbGVhZi4KCmBgYHtyfQpsaWJyYXJ5KGdsbW5ldCkKbGlicmFyeShyZWxhaW1wbykKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKZ2V0IGxlYWZsZW5ndGggZGF0YQpgYGB7cn0KbGVhZmxlbmd0aCA8LSByZWFkX2NzdigiLi4vLi4vcGxhbnQvb3V0cHV0L2xlYWZfbGVuZ3Roc19tZXRhYm9saXRlLmNzdiIpICU+JQogIG11dGF0ZShwb3Q9c3RyX3BhZChwb3QsIHdpZHRoPTMsIHBhZD0iMCIpLAogICAgICAgICBzYW1wbGVJRD1zdHJfYygid3lvIiwgZ2Vub3R5cGUsIHBvdCwgc2VwPSJfIikpICU+JQogIHNlbGVjdChzYW1wbGVJRCwgbGVhZl9hdmdfc3RkKSAKbGVhZmxlbmd0aCAlPiUgYXJyYW5nZShzYW1wbGVJRCkKYGBgCgpnZXQgYW5kIHdyYW5nbGUgbWV0YWJvbGl0ZSBkYXRhCmBgYHtyfQptZXRfcmF3IDwtcmVhZF9jc3YoIi4uL2lucHV0L21ldGFib2xpdGVzX3NldDEuY3N2IikKbWV0IDwtIG1ldF9yYXcgJT4lIAogIG11dGF0ZShwb3Q9c3RyX3BhZChwb3QsIHdpZHRoID0gMywgcGFkID0gIjAiKSkgJT4lCiAgbXV0YXRlKHNhbXBsZUlEPXN0cl9jKCJ3eW8iLCBnZW5vdHlwZSwgcG90LCBzZXA9Il8iKSkgJT4lCiAgc2VsZWN0KHNhbXBsZUlELCBnZW5vdHlwZSwgdGlzc3VlLCBzYW1wbGVfbWFzcyA9IGBzYW1wbGVfbWFzcyBtZ2AsICFzdWJtaXNzaW9uX251bWJlcjpjb25jYXRlbmF0ZSkgJT4lCiAgcGl2b3RfbG9uZ2VyKCFzYW1wbGVJRDpzYW1wbGVfbWFzcywgbmFtZXNfdG8gPSAibWV0YWJvbGl0ZSIsIHZhbHVlc190byA9ICJtZXRfYW1vdW50IikgJT4lCiAgCiAgI2FkanVzdCBieSBzYW1wbGUgbWFzcwogIG11dGF0ZShtZXRfcGVyX21nPW1ldF9hbW91bnQvc2FtcGxlX21hc3MpICU+JQogIAogICNzY2FsZSBhbmQgY2VudGVyCiAgZ3JvdXBfYnkobWV0YWJvbGl0ZSwgZ2Vub3R5cGUsIHRpc3N1ZSkgJT4lCiAgbXV0YXRlKG1ldF9wZXJfbWc9c2NhbGUobWV0X3Blcl9tZyksCiAgICAgICAgIG1ldF9hbXQ9c2NhbGUobWV0X2Ftb3VudCkKICApICU+JSAKICBwaXZvdF93aWRlcihpZF9jb2xzID0gc2FtcGxlSUQsIAogICAgICAgICAgICAgIG5hbWVzX2Zyb20gPSBjKHRpc3N1ZSwgbWV0YWJvbGl0ZSksIAogICAgICAgICAgICAgIHZhbHVlc19mcm9tID0gc3RhcnRzX3dpdGgoIm1ldF8iKSwKICAgICAgICAgICAgICBuYW1lc19zZXAgPSAiXyIpCgptZXQgCmBgYAoKc3BsaXQgdGhpcyBpbnRvIHR3byBkYXRhIGZyYW1lcywgb25lIG5vcm1hbGl6ZWQgYnkgdGlzc3VlIGFtb3VudCBhbmQgb25lIG5vdC4KYGBge3J9Cm1ldF9wZXJfbWcgPC0gbWV0ICU+JSBzZWxlY3Qoc2FtcGxlSUQsICBzdGFydHNfd2l0aCgibWV0X3Blcl9tZyIpKSAlPiUKICBhcy5kYXRhLmZyYW1lKCkgJT4lIGNvbHVtbl90b19yb3duYW1lcygic2FtcGxlSUQiKQptZXRfYW10IDwtIG1ldCAlPiUgc2VsZWN0KHNhbXBsZUlELCAgc3RhcnRzX3dpdGgoIm1ldF9hbXQiKSkgJT4lCiAgYXMuZGF0YS5mcmFtZSgpICU+JSBjb2x1bW5fdG9fcm93bmFtZXMoInNhbXBsZUlEIikKYGBgCgpnZXQgbGVhZiBkYXRhIG9yZGVyIHRvIG1hdGNoCgpgYGB7cn0KbGVhZmxlbmd0aCA8LSBsZWFmbGVuZ3RoW21hdGNoKG1ldCRzYW1wbGVJRCwgbGVhZmxlbmd0aCRzYW1wbGVJRCksXQpsZWFmbGVuZ3RoCmBgYAoKIyMgQ2FsYyBQQ0FzOgoKIyMjIG5vcm1hbGl6ZWQKCiMjIyMgbGVhZgoKYGBge3J9Cm1ldF9wZXJfbWcubGVhZl9QQ0EgPC0gbWV0X3Blcl9tZyAlPiUgCiAgc2VsZWN0KG1hdGNoZXMoIl9sZWFmXyIpKSAlPiUKICBwcmNvbXAoY2VudGVyID0gRkFMU0UsIHNjYWxlLiA9IEZBTFNFKSAjYWxyZWFkeSBjZW50ZXJlZCBhbmQgc2NhbGVkCm5hbWVzKG1ldF9wZXJfbWcubGVhZl9QQ0EpCnRpYmJsZSh2YXJpYW5jZT1tZXRfcGVyX21nLmxlYWZfUENBJHNkZXZeMiwgUEM9c3RyX2MoIlBDIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0cl9wYWQoMTpsZW5ndGgobWV0X3Blcl9tZy5sZWFmX1BDQSRzZGV2KSwgd2lkdGggPSAyLCBwYWQ9IjAiKSkpICU+JQogIG11dGF0ZShwZXJjZW50X3Zhcj0xMDAqdmFyaWFuY2Uvc3VtKHZhcmlhbmNlKSwgIAogICAgICAgICBjdW11bGF0aXZlX3Zhcj1jdW1zdW0ocGVyY2VudF92YXIpKSAlPiUKICBtYWdyaXR0cjo6ZXh0cmFjdCgxOjE1LCkgJT4lCiAgZ2dwbG90KGFlcyh4PVBDLCB5PXBlcmNlbnRfdmFyKSkgKwogIGdlb21fY29sKGZpbGw9InNreWJsdWUiKSArIAogIGdlb21fbGluZShhZXMoeT1jdW11bGF0aXZlX3ZhciksIGdyb3VwPSIiKSArCiAgZ2d0aXRsZSgicGVyY2VudCB2YXJpYW5jZSBleHBsYWluZWQsIG5hbWVkLCBub3JtYWxpemVkIGxlYWYgbWV0YWJvbGl0ZXMiKQpgYGAKCiMjIyMgcm9vdAoKYGBge3J9Cm1ldF9wZXJfbWcucm9vdF9QQ0EgPC0gbWV0X3Blcl9tZyAlPiUgCiAgc2VsZWN0KG1hdGNoZXMoIl9yb290XyIpKSAlPiUKICBwcmNvbXAoY2VudGVyID0gRkFMU0UsIHNjYWxlLiA9IEZBTFNFKSAjYWxyZWFkeSBjZW50ZXJlZCBhbmQgc2NhbGVkCm5hbWVzKG1ldF9wZXJfbWcucm9vdF9QQ0EpCnRpYmJsZSh2YXJpYW5jZT1tZXRfcGVyX21nLnJvb3RfUENBJHNkZXZeMiwgUEM9c3RyX2MoIlBDIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0cl9wYWQoMTpsZW5ndGgobWV0X3Blcl9tZy5yb290X1BDQSRzZGV2KSwgd2lkdGggPSAyLCBwYWQ9IjAiKSkpICU+JQogIG11dGF0ZShwZXJjZW50X3Zhcj0xMDAqdmFyaWFuY2Uvc3VtKHZhcmlhbmNlKSwgIAogICAgICAgICBjdW11bGF0aXZlX3Zhcj1jdW1zdW0ocGVyY2VudF92YXIpKSAlPiUKICBtYWdyaXR0cjo6ZXh0cmFjdCgxOjE1LCkgJT4lCiAgZ2dwbG90KGFlcyh4PVBDLCB5PXBlcmNlbnRfdmFyKSkgKwogIGdlb21fY29sKGZpbGw9InNreWJsdWUiKSArIAogIGdlb21fbGluZShhZXMoeT1jdW11bGF0aXZlX3ZhciksIGdyb3VwPSIiKSArCiAgZ2d0aXRsZSgicGVyY2VudCB2YXJpYW5jZSBleHBsYWluZWQsIG5hbWVkLCBub3JtYWxpemVkIHJvb3QgbWV0YWJvbGl0ZXMiKQpgYGAKCiMjIyByYXcKCiMjIyMgbGVhZgpgYGB7cn0KbWV0X2FtdC5sZWFmX1BDQSA8LSBtZXRfYW10ICU+JQogIHNlbGVjdChtYXRjaGVzKCJfbGVhZl8iKSkgJT4lCiAgcHJjb21wKGNlbnRlciA9IEZBTFNFLCBzY2FsZS4gPSBGQUxTRSkgI2FscmVhZHkgY2VudGVyZWQgYW5kIHNjYWxlZApuYW1lcyhtZXRfcGVyX21nLmxlYWZfUENBKQp0aWJibGUodmFyaWFuY2U9bWV0X2FtdC5sZWFmX1BDQSRzZGV2XjIsIFBDPXN0cl9jKCJQQyIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdHJfcGFkKDE6bGVuZ3RoKG1ldF9hbXQubGVhZl9QQ0Ekc2RldiksIHdpZHRoID0gMiwgcGFkPSIwIikpKSAlPiUKICBtdXRhdGUocGVyY2VudF92YXI9MTAwKnZhcmlhbmNlL3N1bSh2YXJpYW5jZSksICAKICAgICAgICAgY3VtdWxhdGl2ZV92YXI9Y3Vtc3VtKHBlcmNlbnRfdmFyKSkgJT4lCiAgbWFncml0dHI6OmV4dHJhY3QoMToxNSwpICU+JQogIGdncGxvdChhZXMoeD1QQywgeT1wZXJjZW50X3ZhcikpICsKICBnZW9tX2NvbChmaWxsPSJza3libHVlIikgKyAKICBnZW9tX2xpbmUoYWVzKHk9Y3VtdWxhdGl2ZV92YXIpLCBncm91cD0iIikgKwogIGdndGl0bGUoInBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkLCBuYW1lZCwgcmF3IGxlYWYgbWV0YWJvbGl0ZXMiKQpgYGAKCiMjIyMgcm9vdApgYGB7cn0KbWV0X2FtdC5yb290X1BDQSA8LSBtZXRfYW10ICU+JQogIHNlbGVjdChtYXRjaGVzKCJfcm9vdF8iKSkgJT4lCiAgcHJjb21wKGNlbnRlciA9IEZBTFNFLCBzY2FsZS4gPSBGQUxTRSkgI2FscmVhZHkgY2VudGVyZWQgYW5kIHNjYWxlZApuYW1lcyhtZXRfcGVyX21nLnJvb3RfUENBKQp0aWJibGUodmFyaWFuY2U9bWV0X2FtdC5yb290X1BDQSRzZGV2XjIsIFBDPXN0cl9jKCJQQyIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdHJfcGFkKDE6bGVuZ3RoKG1ldF9hbXQucm9vdF9QQ0Ekc2RldiksIHdpZHRoID0gMiwgcGFkPSIwIikpKSAlPiUKICBtdXRhdGUocGVyY2VudF92YXI9MTAwKnZhcmlhbmNlL3N1bSh2YXJpYW5jZSksICAKICAgICAgICAgY3VtdWxhdGl2ZV92YXI9Y3Vtc3VtKHBlcmNlbnRfdmFyKSkgJT4lCiAgbWFncml0dHI6OmV4dHJhY3QoMToxNSwpICU+JQogIGdncGxvdChhZXMoeD1QQywgeT1wZXJjZW50X3ZhcikpICsKICBnZW9tX2NvbChmaWxsPSJza3libHVlIikgKyAKICBnZW9tX2xpbmUoYWVzKHk9Y3VtdWxhdGl2ZV92YXIpLCBncm91cD0iIikgKwogIGdndGl0bGUoInBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkLCBuYW1lZCwgcmF3IHJvb3QgbWV0YWJvbGl0ZXMiKQpgYGAKCiMjIG5vdyB0cnkgdGhlc2UgaW4gYSBwZW5hbGl6ZWQgcmVncmVzc2lvbgoKIyMjIG5vcm1hbGl6ZWQKCmFyZSB0aGUgUENzIG5vcm1hbGl6ZWQ/CmBgYHtyfQpjb2xNZWFucyhtZXRfYW10LmxlYWZfUENBJHgpICU+JSByb3VuZCgzKSAjeWVzIGNlbnRlcmVkCmFwcGx5KG1ldF9hbXQubGVhZl9QQ0EkeCwgMiwgc2QpICU+JSByb3VuZCgyKSAjbm90IHNjYWxlZApgYGAKCmNvbWJpbmUgdGhlIGxlYWYgYW5kIHJvb3QsIGFuZCB0aGVuIHNjYWxlIHRoZW06CmBgYHtyfQptZXRfcGVyX21nLmxlYWZfUENzIDwtIG1ldF9wZXJfbWcubGVhZl9QQ0EkeApjb2xuYW1lcyhtZXRfcGVyX21nLmxlYWZfUENzKSA8LSBzdHJfYygibGVhZl8iLCBjb2xuYW1lcyhtZXRfcGVyX21nLmxlYWZfUENzKSkKCm1ldF9wZXJfbWcucm9vdF9QQ3MgPC0gbWV0X3Blcl9tZy5yb290X1BDQSR4CmNvbG5hbWVzKG1ldF9wZXJfbWcucm9vdF9QQ3MpIDwtIHN0cl9jKCJyb290XyIsIGNvbG5hbWVzKG1ldF9wZXJfbWcucm9vdF9QQ3MpKQoKbWV0X3Blcl9tZy5QQ3MgPC0gY2JpbmQobWV0X3Blcl9tZy5sZWFmX1BDcywgbWV0X3Blcl9tZy5yb290X1BDcykgJT4lCiAgc2NhbGUoKQoKbWV0X2FtdC5sZWFmX1BDcyA8LSBtZXRfYW10LmxlYWZfUENBJHgKY29sbmFtZXMobWV0X2FtdC5sZWFmX1BDcykgPC0gc3RyX2MoImxlYWZfIiwgY29sbmFtZXMobWV0X2FtdC5sZWFmX1BDcykpCgptZXRfYW10LnJvb3RfUENzIDwtIG1ldF9hbXQucm9vdF9QQ0EkeApjb2xuYW1lcyhtZXRfYW10LnJvb3RfUENzKSA8LSBzdHJfYygicm9vdF8iLCBjb2xuYW1lcyhtZXRfYW10LnJvb3RfUENzKSkKCm1ldF9hbXQuUENzIDwtIGNiaW5kKG1ldF9hbXQubGVhZl9QQ3MsIG1ldF9hbXQucm9vdF9QQ3MpICU+JQogIHNjYWxlKCkKYGBgCgphbHNvIGNvbWJpbmUgdGhlIHJvdGF0aW9ucwpgYGB7cn0KbWV0X3Blcl9tZy5sZWFmX3JvdGF0aW9uIDwtIG1ldF9wZXJfbWcubGVhZl9QQ0Ekcm90YXRpb24gJT4lCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICByZW5hbWVfd2l0aCh+IHN0cl9jKCJsZWFmXyIsIC54KSkgJT4lCiAgcm93bmFtZXNfdG9fY29sdW1uKCJtZXRhYm9saXRlIikKCm1ldF9wZXJfbWcucm9vdF9yb3RhdGlvbiA8LSBtZXRfcGVyX21nLnJvb3RfUENBJHJvdGF0aW9uICU+JQogIGFzLmRhdGEuZnJhbWUoKSAlPiUgCiAgcmVuYW1lX3dpdGgofiBzdHJfYygicm9vdF8iLCAueCkpICU+JQogIHJvd25hbWVzX3RvX2NvbHVtbigibWV0YWJvbGl0ZSIpCgptZXRfcGVyX21nLlBDX3JvdGF0aW9uIDwtIGZ1bGxfam9pbihtZXRfcGVyX21nLmxlYWZfcm90YXRpb24sIG1ldF9wZXJfbWcucm9vdF9yb3RhdGlvbiwgYnk9Im1ldGFib2xpdGUiKQoKbWV0X2FtdC5sZWFmX3JvdGF0aW9uIDwtIG1ldF9hbXQubGVhZl9QQ0Ekcm90YXRpb24gJT4lIAogIGFzLmRhdGEuZnJhbWUoKSAlPiUgCiAgcmVuYW1lX3dpdGgofiBzdHJfYygibGVhZl8iLCAueCkpICU+JQogIHJvd25hbWVzX3RvX2NvbHVtbigibWV0YWJvbGl0ZSIpCgptZXRfYW10LnJvb3Rfcm90YXRpb24gPC0gbWV0X2FtdC5yb290X1BDQSRyb3RhdGlvbiAlPiUKICBhcy5kYXRhLmZyYW1lKCkgJT4lIAogIHJlbmFtZV93aXRoKH4gc3RyX2MoInJvb3RfIiwgLngpKSAlPiUKICByb3duYW1lc190b19jb2x1bW4oIm1ldGFib2xpdGUiKQoKbWV0X2FtdC5QQ19yb3RhdGlvbiA8LSBmdWxsX2pvaW4obWV0X2FtdC5sZWFmX3JvdGF0aW9uLCBtZXRfYW10LnJvb3Rfcm90YXRpb24sIGJ5PSJtZXRhYm9saXRlIikKCmBgYAoKCgpgYGB7cn0KbWV0X3Blcl9tZ19maXQxTE9PIDwtIGN2LmdsbW5ldCh4PW1ldF9wZXJfbWcuUENzLCB5PWxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkLCBuZm9sZHMgPSBucm93KG1ldF9wZXJfbWcuUENzKSwgYWxwaGE9MSApCnBsb3QobWV0X3Blcl9tZ19maXQxTE9PKQpiZXN0bGFtPW1ldF9wZXJfbWdfZml0MUxPTyRsYW1iZGEuMXNlCmBgYAoKICAKCk5FWFQgU1RFUDogRG8gYSBLLWZvbGQgQ1YsIHJlcGVhdCBtYW55IHRpbWVzIGFuZCBhdmVyYWdlLiAgTWlnaHQgYXMgd2VsbCBkbyBhbHBoYSB3aGlsZSB3ZSBhcmUgYXQgaXQuICBJZiB3ZSBhcmUgZG9pbmcgYWxwaGEsIHRoZW4gd2UgbmVlZCB0byBtYW51YWxseSBjcmVhdGUgb3VyIG93biBmb2xkcyBsaXN0IGZvciBlYWNoIHJ1bgoKIyBub3JtYWxpemVkCgojIyBtdWx0aSBDVgoKRml0IDEwMSBDVnMgZm9yIGVhY2ggb2YgMTEgYWxwaGFzCmBgYHtyfQpzZXQuc2VlZCgxMjQ1KQoKZm9sZHMgPC0gdGliYmxlKHJ1bj0xOjEwMSkgJT4lIAogIG11dGF0ZShmb2xkcz1tYXAocnVuLCB+IHNhbXBsZShyZXAoMTo2LDYpKSkpCgpzeXN0ZW0udGltZSAobWV0X3Blcl9tZ19tdWx0aUNWIDwtIGV4cGFuZF9ncmlkKHJ1bj0xOjEwMCwgYWxwaGE9cm91bmQoc2VxKDAsMSwuMSksMSkpICU+JQogICAgICAgICAgICAgICBsZWZ0X2pvaW4oZm9sZHMsIGJ5PSJydW4iKSAlPiUKICAgICAgICAgICAgICAgbXV0YXRlKGZpdD1tYXAyKGZvbGRzLCBhbHBoYSwgfiBjdi5nbG1uZXQoeD1tZXRfcGVyX21nLlBDcywgeT1sZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCwgZm9sZGlkID0gLngsIGFscGhhPS55CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICkpKQogICAgICAgICAgICAgIywgbGFtYmRhPWV4cChzZXEoLTUsMCxsZW5ndGgub3V0ID0gNTApKSApKSkKKSAjMTAwIHNlY29uZHMKCmhlYWQobWV0X3Blcl9tZ19tdWx0aUNWKQpgYGAKCmZvciBlYWNoIGZpdCwgcHVsbCBvdXQgdGhlIG1lYW4gY3YgZXJyb3IsIGxhbWJkYSwgbWluIGxhbWJkYSwgYW5kIDFzZSBsYW1iZGEgCmBgYHtyfQptZXRfcGVyX21nX211bHRpQ1YgPC0gbWV0X3Blcl9tZ19tdWx0aUNWICU+JQogIG11dGF0ZShjdm09bWFwKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoImN2bSIpKSwKICAgICAgICAgbGFtYmRhPW1hcChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJsYW1iZGEiKSksCiAgICAgICAgIGxhbWJkYS5taW49bWFwX2RibChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJsYW1iZGEubWluIiApKSwKICAgICAgICAgbGFtYmRhLjFzZT1tYXBfZGJsKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoImxhbWJkYS4xc2UiKSksCiAgICAgICAgIG56ZXJvPW1hcChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJuemVybyIpKQogICkKCmhlYWQobWV0X3Blcl9tZ19tdWx0aUNWKQpgYGAKCgpub3cgY2FsY3VsYXRlIHRoZSBtZWFuIGFuZCBzZW0gb2YgY3ZtIGFuZCBtaW4sMXNlIGxhYm1kYXMuICBUaGVzZSBuZWVkIHRvIGJlIGRvbmUgc2VwYXJhdGVseSBiZWNhdXNlIG9mIHRoZSB3YXkgdGhlIGdyb3VwaW5nIHdvcmtzCmBgYHtyfQptZXRfcGVyX21nX3N1bW1hcnlfY3ZtIDwtIG1ldF9wZXJfbWdfbXVsdGlDViAlPiUgZHBseXI6OnNlbGVjdCgtZml0LCAtZm9sZHMpICU+JSAKICB1bm5lc3QoYyhjdm0sIGxhbWJkYSkpICU+JQogIGdyb3VwX2J5KGFscGhhLCBsYW1iZGEpICU+JQogIHN1bW1hcml6ZShtZWFuY3ZtPW1lYW4oY3ZtKSwgc2VtPXNkKGN2bSkvc3FydChuKCkpLCBoaWdoPW1lYW5jdm0rc2VtLCBsb3c9bWVhbmN2bS1zZW0pCgptZXRfcGVyX21nX3N1bW1hcnlfY3ZtCmBgYAoKYGBge3J9Cm1ldF9wZXJfbWdfc3VtbWFyeV9sYW1iZGEgPC0gbWV0X3Blcl9tZ19tdWx0aUNWICU+JSBkcGx5cjo6c2VsZWN0KC1maXQsIC1mb2xkcywgLWN2bSkgJT4lIAogIGdyb3VwX2J5KGFscGhhKSAlPiUKICBzdW1tYXJpemUoCiAgICBsYW1iZGEubWluLnNkPXNkKGxhbWJkYS5taW4pLCAKICAgIGxhbWJkYS5taW4ubWVhbj1tZWFuKGxhbWJkYS5taW4pLAogICAgI2xhbWJkYS5taW4ubWVkPW1lZGlhbihsYW1iZGEubWluKSwgCiAgICBsYW1iZGEubWluLmhpZ2g9bGFtYmRhLm1pbi5tZWFuK2xhbWJkYS5taW4uc2QsCiAgICAjbGFtYmRhLm1pbi5sb3c9bGFtYmRhLm1pbi5tZWFuLWxhbWJkYS5taW4uc2VtLAogICAgI2xhbWJkYS4xc2Uuc2VtPXNkKGxhbWJkYS4xc2UpL3NxcnQobigpKSwgCiAgICBsYW1iZGEuMXNlLm1lYW49bWVhbihsYW1iZGEuMXNlKSwKICAgICNsYW1iZGEuMXNlLm1lZD1tZWRpYW4obGFtYmRhLjFzZSksIAogICAgI2xhbWJkYS4xc2UuaGlnaD1sYW1iZGEuMXNlK2xhbWJkYS4xc2Uuc2VtLAogICAgI2xhbWJkYS4xc2UubG93PWxhbWJkYS4xc2UtbGFtYmRhLjFzZS5zZW0sCiAgICBuemVybz1uemVyb1sxXSwKICAgIGxhbWJkYT1sYW1iZGFbMV0KICApCgptZXRfcGVyX21nX3N1bW1hcnlfbGFtYmRhCmBgYAoKCnBsb3QgaXQKYGBge3J9Cm1ldF9wZXJfbWdfc3VtbWFyeV9jdm0gJT4lCiAgI2ZpbHRlcihhbHBoYSE9MCkgJT4lICMgd29yc2UgdGhhbiBldmVyeXRoaW5nIGVsc2UgYW5kIHRocm93aW5nIHRoZSBwbG90cyBvZmYKICBnZ3Bsb3QoYWVzKHg9bG9nKGxhbWJkYSksIHk9IG1lYW5jdm0sICB5bWluPWxvdywgeW1heD1oaWdoKSkgKwogIGdlb21fcmliYm9uKGFscGhhPS4yNSkgKwogIGdlb21fbGluZShhZXMoY29sb3I9YXMuY2hhcmFjdGVyKGFscGhhKSkpICsKICBmYWNldF93cmFwKH4gYXMuY2hhcmFjdGVyKGFscGhhKSkgKwogICBjb29yZF9jYXJ0ZXNpYW4oeGxpbT0oYygtNSwwKSkpICsKICBnZW9tX3ZsaW5lKGFlcyh4aW50ZXJjZXB0PWxvZyhsYW1iZGEubWluLm1lYW4pKSwgYWxwaGE9LjUsIGRhdGE9bWV0X3Blcl9tZ19zdW1tYXJ5X2xhbWJkYSkgKwogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQ9bG9nKGxhbWJkYS5taW4uaGlnaCkpLCBhbHBoYT0uNSwgZGF0YT1tZXRfcGVyX21nX3N1bW1hcnlfbGFtYmRhLCBjb2xvcj0iYmx1ZSIpIAoKYGBgCgoKClNvIG92ZXJhbGwgdGhlc2UgbG9vayBtb3JlIHJlYXNvbmFibGUgdGhhbiB0aGUgTE9PIHBsb3QuCgpNYWtlIGEgcGxvdCBvZiBNU0UgYXQgbWluaW11bSBsYW1iZGEgZm9yIGVhY2ggYWxwaGEKCmBgYHtyfQptZXRfcGVyX21nX3N1bW1hcnlfY3ZtICU+JSAKICBncm91cF9ieShhbHBoYSkgJT4lCiAgZmlsdGVyKHJhbmsobWVhbmN2bSwgdGllcy5tZXRob2QgPSAiZmlyc3QiKT09MSkgJT4lCiAgZ2dwbG90KGFlcyh4PWFscGhhLHk9bWVhbmN2bSx5bWluPWxvdyx5bWF4PWhpZ2gpKSArCiAgZ2VvbV9yaWJib24oY29sb3I9TkEsIGZpbGw9ImdyYXk4MCIpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludCgpCmBgYApub3QgYSBwYXJ0aWN1bGFyIGxhcmdlIGRpZmZlcmVuY2UgdGhlcmUsIGFzaWRlIGZyb20gMC4xIGFuZCBldmVuIHRoZW4sIG5vdCB0b28gbXVjaCBiZXR0ZXIuCgpQbG90IHRoZSBudW1iZXIgb2Ygbnplcm8gY29lZmZpY2llbnRzCgpgYGB7cn0KbWV0X3Blcl9tZ19zdW1tYXJ5X2xhbWJkYSAlPiUKICB1bm5lc3QoYyhsYW1iZGEsIG56ZXJvKSkgJT4lCiAgZ3JvdXBfYnkoYWxwaGEpICU+JQogIGZpbHRlcihhYnMobGFtYmRhLm1pbi5tZWFuLWxhbWJkYSk9PW1pbihhYnMobGFtYmRhLm1pbi5tZWFuLWxhbWJkYSkpICApICU+JQogIHVuZ3JvdXAoKSAlPiUKCmdncGxvdChhZXMoeD1hcy5jaGFyYWN0ZXIoYWxwaGEpLCB5PW56ZXJvKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2d0aXRsZSgiTnVtYmVyIG9mIG5vbi16ZXJvIGNvZWZmaWNlbnRzIGF0IG1pbmltdW0gbGFtYmRhIikgKwogIHlsaW0oMCwzNikKYGBgCk9LIGxldCdzIGRvIHJlcGVhdGVkIHRlc3QgdHJhaW4gc3RhcnRpbmcgZnJvbSB0aGVzZSBDViBsYW1iZGFzCgpgYGB7cn0KbXVsdGlfdHQgPC0gZnVuY3Rpb24obGFtYmRhLCBhbHBoYSwgbj0xMDAwMCwgc2FtcGxlX3NpemU9MzYsIHRyYWluX3NpemU9MzAsIHgsIHk9bGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQpIHsKICBwcmludChsYW1iZGEpCiAgcHJpbnQoYWxwaGEpCnR0IDwtCiAgdGliYmxlKHJ1bj0xOm4pICU+JQogIG11dGF0ZSh0cmFpbj1tYXAocnVuLCB+IHNhbXBsZSgxOnNhbXBsZV9zaXplLCB0cmFpbl9zaXplKSkpICU+JQogIG11dGF0ZShmaXQ9bWFwKHRyYWluLCB+IGdsbW5ldCh4PXhbLixdLCB5PXlbLl0sIGxhbWJkYSA9IGxhbWJkYSwgYWxwaGEgPSBhbHBoYSApKSkgJT4lCiAgCiAgbXV0YXRlKHByZWQ9bWFwMihmaXQsIHRyYWluLCB+IHByZWRpY3QoLngsIG5ld3ggPSB4Wy0ueSxdKSkpICU+JQogIG11dGF0ZShjb3I9bWFwMl9kYmwocHJlZCwgdHJhaW4sIH4gY29yKC54LCB5Wy0ueV0pICApKSAlPiUKICBtdXRhdGUoTVNFPW1hcDJfZGJsKHByZWQsIHRyYWluLCB+IG1lYW4oKHlbLS55XSAtIC54KV4yKSkpICU+JQogIHN1bW1hcml6ZSgKICAgIG51bV9uYT1zdW0oaXMubmEoY29yKSksIAogICAgbnVtX2x0XzA9c3VtKGNvcjw9MCwgbmEucm09VFJVRSksCiAgICBhdmdfY29yPW1lYW4oY29yLCBuYS5ybT1UUlVFKSwKICAgIGF2Z19NU0U9bWVhbihNU0UpKQp0dAp9CgpwZXJfbWdfZml0X3Rlc3RfdHJhaW4gPC0gbWV0X3Blcl9tZ19zdW1tYXJ5X2xhbWJkYSAlPiUgCiAgc2VsZWN0KGFscGhhLCBsYW1iZGEubWluLm1lYW4pCgpwZXJfbWdfZml0X3Rlc3RfdHJhaW4gPC0gbWV0X3Blcl9tZ19tdWx0aUNWICU+JQogIGZpbHRlcihydW49PTEpICU+JQogIHNlbGVjdChhbHBoYSwgZml0KSAlPiUKICByaWdodF9qb2luKHBlcl9tZ19maXRfdGVzdF90cmFpbikKCnBlcl9tZ19maXRfdGVzdF90cmFpbiA8LSBwZXJfbWdfZml0X3Rlc3RfdHJhaW4gJT4lCiAgbXV0YXRlKHByZWRfZnVsbD1tYXAyKGZpdCwgbGFtYmRhLm1pbi5tZWFuLCB+IHByZWRpY3QoLngsIHM9LnksIG5ld3g9bWV0X3Blcl9tZy5QQ3MpKSwKICAgICAgICAgZnVsbF9SPW1hcF9kYmwocHJlZF9mdWxsLCB+IGNvcigueCwgbGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQpKSwKICAgICAgICAgZnVsbF9NU0U9bWFwX2RibChwcmVkX2Z1bGwsIH4gbWVhbigobGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQtLngpXjIpKSkgJT4lCiAgCiAgbXV0YXRlKHR0PW1hcDIobGFtYmRhLm1pbi5tZWFuLCBhbHBoYSwgfiBtdWx0aV90dChsYW1iZGE9LngsIGFscGhhPS55LCB4PW1ldF9wZXJfbWcuUENzKSkpCgoKCihwZXJfbWdfZml0X3Rlc3RfdHJhaW4gPC0gcGVyX21nX2ZpdF90ZXN0X3RyYWluICU+JSB1bm5lc3QodHQpKQpgYGAKCmBgYHtyfQpwZXJfbWdfZml0X3Rlc3RfdHJhaW4gJT4lCiAgZ2dwbG90KGFlcyh4PWFscGhhKSkgKwogIGdlb21fbGluZShhZXMoeT1hdmdfY29yKSwgY29sb3I9InJlZCIpICsKICBnZW9tX3BvaW50KGFlcyh5PWF2Z19jb3IpLCBjb2xvcj0icmVkIikgKwogIGdlb21fbGluZShhZXMoeT1hdmdfTVNFKSwgY29sb3I9ImJsdWUiKSArCiAgZ2VvbV9wb2ludChhZXMoeT1hdmdfTVNFKSwgY29sb3I9ImJsdWUiKQpgYGAKYWxwaGEgb2YgMC44IHRvIDEuMCBhcmUgdmVyeSBzaW1pbGFyIGFuZCBhcmUgdGhlIGJlc3QgaGVyZS4KCiMjIGxvb2sgYXQgZml0OgoKYGBge3J9CmFscGhhX3Blcl9tZyA8LSAuOAoKYmVzdF9wZXJfbWcgPC0gcGVyX21nX2ZpdF90ZXN0X3RyYWluICU+JSBmaWx0ZXIoYWxwaGEgPT0gYWxwaGFfcGVyX21nKSAKYmVzdF9wZXJfbWdfZml0IDwtIGJlc3RfcGVyX21nJGZpdFtbMV1dCmJlc3RfcGVyX21nX2xhbWJkYSA8LSBiZXN0X3Blcl9tZyRsYW1iZGEubWluLm1lYW4KCnBlcl9tZ19jb2VmLnRiIDwtIGNvZWYoYmVzdF9wZXJfbWdfZml0LCBzPWJlc3RfcGVyX21nX2xhbWJkYSkgJT4lIAogIGFzLm1hdHJpeCgpICU+JSBhcy5kYXRhLmZyYW1lKCkgJT4lIAogIHJvd25hbWVzX3RvX2NvbHVtbih2YXI9IlBDIikgJT4lCiAgcmVuYW1lKGJldGE9YDFgKQogIApwZXJfbWdfY29lZi50YiAlPiUgZmlsdGVyKGJldGEhPTApICU+JSBhcnJhbmdlKGJldGEpCgpgYGAKCnByZWQgYW5kIG9icwpgYGB7cn0KcGxvdChsZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCwgYmVzdF9wZXJfbWckcHJlZF9mdWxsW1sxXV0pCmNvci50ZXN0KGxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkLCBiZXN0X3Blcl9tZyRwcmVkX2Z1bGxbWzFdXSkgIy41NwpiZXN0X3Blcl9tZyRmdWxsX01TRQpgYGAKCiMjIFBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkCgpgYGB7cn0KcGVyX21nX3ZhcnMgPC0gcGVyX21nX2NvZWYudGIgJT4lIAogIGZpbHRlcihiZXRhICE9MCwgUEMhPSIoSW50ZXJjZXB0KSIpICU+JQogIHB1bGwoUEMpICU+JSBjKCJsZWFmX2F2Z19zdGQiLCAuKQoKcGVyX21nX3JlbGltcCA8LSBsZWFmbGVuZ3RoICU+JSBzZWxlY3QobGVhZl9hdmdfc3RkKSAlPiUgY2JpbmQobWV0X3Blcl9tZy5QQ3MpICU+JSBhcy5kYXRhLmZyYW1lKCkgJT4lIGRwbHlyOjpzZWxlY3QoYWxsX29mKHBlcl9tZ192YXJzKSkgJT4lCiAgY2FsYy5yZWxpbXAoKSAKCnBlcl9tZ19jb2VmLnRiIDwtIHBlcl9tZ19yZWxpbXBAbG1nICU+JSBhcy5tYXRyaXgoKSAlPiUgYXMuZGF0YS5mcmFtZSgpICU+JQogIHJvd25hbWVzX3RvX2NvbHVtbigiUEMiKSAlPiUKICByZW5hbWUoUHJvcFZhcl9tZXRfcGVyX21nPVYxKSAlPiUKICBmdWxsX2pvaW4ocGVyX21nX2NvZWYudGIpICU+JQogIGFycmFuZ2UoZGVzYyhQcm9wVmFyX21ldF9wZXJfbWcpKQoKcGVyX21nX2NvZWYudGIKCmBgYAoKQ2hlY2tvdXQgdGhlIHJvdGF0aW9ucy4gIAoKCmBgYHtyfQptZXRfcGVyX21nX3JvdGF0aW9uX291dCA8LSBtZXRfcGVyX21nLlBDX3JvdGF0aW9uICU+JSAKICBwaXZvdF9sb25nZXIoLW1ldGFib2xpdGUsIG5hbWVzX3RvPSJQQyIsIHZhbHVlc190bz0ibG9hZGluZyIpICU+JQogIGZpbHRlcihQQyAlaW4lIGZpbHRlcihwZXJfbWdfY29lZi50YiwgYmV0YSE9MCkkUEMgKSAlPiUKICBncm91cF9ieShQQykgJT4lCiAgICBmaWx0ZXIoIXN0cl9kZXRlY3QobWV0YWJvbGl0ZSwiLioobGVhZnxyb290KV9bMC05XSokIikpICU+JQogIGZpbHRlcihhYnMobG9hZGluZykgPj0gMC4wNSkgJT4lCiAgbGVmdF9qb2luKHBlcl9tZ19jb2VmLnRiLCBieT0iUEMiKSAlPiUKICBhcnJhbmdlKGRlc2MoYWJzKGJldGEpKSwgZGVzYyhhYnMobG9hZGluZykpKSAlPiUKICBtdXRhdGUob3JnYW49aWZlbHNlKHN0cl9kZXRlY3QobWV0YWJvbGl0ZSwgIl9sZWFmXyIpLCAibGVhZiIsICJyb290IiksCiAgICAgICAgIHRyYW5zZm9ybWF0aW9uPSJub3JtYWxpemVkIiwKICAgICAgICAgbWV0YWJvbGl0ZT1zdHJfcmVtb3ZlKG1ldGFib2xpdGUsICJtZXRfcGVyX21nXyhyb290fGxlYWYpXyIpLAogICAgICAgICBtZXRhYm9saXRlX2VmZmVjdF9vbl9sZWFmPWlmZWxzZShiZXRhKmxvYWRpbmc+MCwgImluY3JlYXNlIiwgImRlY3JlYXNlIikpCm1ldF9wZXJfbWdfcm90YXRpb25fb3V0ICU+JSAgd3JpdGVfY3N2KCIuLi9vdXRwdXQvTGVhZl9hc3NvY2lhdGVkX21ldGFib2xpdGVzX25vcm1hbGl6ZWQuY3N2IikKCm1ldF9wZXJfbWdfcm90YXRpb25fb3V0CmBgYAoKIyBub24tbm9ybWF6bGl6ZWQKCiMjIG11bHRpIENWCgpGaXQgMTAxIENWcyBmb3IgZWFjaCBvZiAxMSBhbHBoYXMKYGBge3J9CnNldC5zZWVkKDEyNDUpCgpmb2xkcyA8LSB0aWJibGUocnVuPTE6MTAxKSAlPiUgCiAgbXV0YXRlKGZvbGRzPW1hcChydW4sIH4gc2FtcGxlKHJlcCgxOjYsNikpKSkKCnN5c3RlbS50aW1lIChtZXRfYW10X211bHRpQ1YgPC0gZXhwYW5kX2dyaWQocnVuPTE6MTAwLCBhbHBoYT1yb3VuZChzZXEoMCwxLC4xKSwxKSkgJT4lCiAgICAgICAgICAgICAgIGxlZnRfam9pbihmb2xkcywgYnk9InJ1biIpICU+JQogICAgICAgICAgICAgICBtdXRhdGUoZml0PW1hcDIoZm9sZHMsIGFscGhhLCB+IGN2LmdsbW5ldCh4PW1ldF9hbXQuUENzLCB5PWxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkLCBmb2xkaWQgPSAueCwgYWxwaGE9LnkKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgKSkpCiAgICAgICAgICAgICAjLCBsYW1iZGE9ZXhwKHNlcSgtNSwwLGxlbmd0aC5vdXQgPSA1MCkpICkpKQopICMxMDAgc2Vjb25kcwoKaGVhZChtZXRfYW10X211bHRpQ1YpCmBgYAoKZm9yIGVhY2ggZml0LCBwdWxsIG91dCB0aGUgbWVhbiBjdiBlcnJvciwgbGFtYmRhLCBtaW4gbGFtYmRhLCBhbmQgMXNlIGxhbWJkYSAKYGBge3J9Cm1ldF9hbXRfbXVsdGlDViA8LSBtZXRfYW10X211bHRpQ1YgJT4lCiAgbXV0YXRlKGN2bT1tYXAoZml0LCBtYWdyaXR0cjo6ZXh0cmFjdCgiY3ZtIikpLAogICAgICAgICBsYW1iZGE9bWFwKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoImxhbWJkYSIpKSwKICAgICAgICAgbGFtYmRhLm1pbj1tYXBfZGJsKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoImxhbWJkYS5taW4iICkpLAogICAgICAgICBsYW1iZGEuMXNlPW1hcF9kYmwoZml0LCBtYWdyaXR0cjo6ZXh0cmFjdCgibGFtYmRhLjFzZSIpKSwKICAgICAgICAgbnplcm89bWFwKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoIm56ZXJvIikpCiAgKQoKaGVhZChtZXRfYW10X211bHRpQ1YpCmBgYAoKCm5vdyBjYWxjdWxhdGUgdGhlIG1lYW4gYW5kIHNlbSBvZiBjdm0gYW5kIG1pbiwxc2UgbGFibWRhcy4gIFRoZXNlIG5lZWQgdG8gYmUgZG9uZSBzZXBhcmF0ZWx5IGJlY2F1c2Ugb2YgdGhlIHdheSB0aGUgZ3JvdXBpbmcgd29ya3MKYGBge3J9Cm1ldF9hbXRfc3VtbWFyeV9jdm0gPC0gbWV0X2FtdF9tdWx0aUNWICU+JSBkcGx5cjo6c2VsZWN0KC1maXQsIC1mb2xkcykgJT4lIAogIHVubmVzdChjKGN2bSwgbGFtYmRhKSkgJT4lCiAgZ3JvdXBfYnkoYWxwaGEsIGxhbWJkYSkgJT4lCiAgc3VtbWFyaXplKG1lYW5jdm09bWVhbihjdm0pLCBzZW09c2QoY3ZtKS9zcXJ0KG4oKSksIGhpZ2g9bWVhbmN2bStzZW0sIGxvdz1tZWFuY3ZtLXNlbSkKCm1ldF9hbXRfc3VtbWFyeV9jdm0KYGBgCgpgYGB7cn0KbWV0X2FtdF9zdW1tYXJ5X2xhbWJkYSA8LSBtZXRfYW10X211bHRpQ1YgJT4lIGRwbHlyOjpzZWxlY3QoLWZpdCwgLWZvbGRzLCAtY3ZtKSAlPiUgCiAgZ3JvdXBfYnkoYWxwaGEpICU+JQogIHN1bW1hcml6ZSgKICAgIGxhbWJkYS5taW4uc2Q9c2QobGFtYmRhLm1pbiksIAogICAgbGFtYmRhLm1pbi5tZWFuPW1lYW4obGFtYmRhLm1pbiksCiAgICAjbGFtYmRhLm1pbi5tZWQ9bWVkaWFuKGxhbWJkYS5taW4pLCAKICAgIGxhbWJkYS5taW4uaGlnaD1sYW1iZGEubWluLm1lYW4rbGFtYmRhLm1pbi5zZCwKICAgICNsYW1iZGEubWluLmxvdz1sYW1iZGEubWluLm1lYW4tbGFtYmRhLm1pbi5zZW0sCiAgICAjbGFtYmRhLjFzZS5zZW09c2QobGFtYmRhLjFzZSkvc3FydChuKCkpLCAKICAgIGxhbWJkYS4xc2UubWVhbj1tZWFuKGxhbWJkYS4xc2UpLAogICAgI2xhbWJkYS4xc2UubWVkPW1lZGlhbihsYW1iZGEuMXNlKSwgCiAgICAjbGFtYmRhLjFzZS5oaWdoPWxhbWJkYS4xc2UrbGFtYmRhLjFzZS5zZW0sCiAgICAjbGFtYmRhLjFzZS5sb3c9bGFtYmRhLjFzZS1sYW1iZGEuMXNlLnNlbSwKICAgIG56ZXJvPW56ZXJvWzFdLAogICAgbGFtYmRhPWxhbWJkYVsxXQogICkKCm1ldF9hbXRfc3VtbWFyeV9sYW1iZGEKYGBgCgoKcGxvdCBpdApgYGB7cn0KbWV0X2FtdF9zdW1tYXJ5X2N2bSAlPiUKICAjZmlsdGVyKGFscGhhIT0wKSAlPiUgIyB3b3JzZSB0aGFuIGV2ZXJ5dGhpbmcgZWxzZSBhbmQgdGhyb3dpbmcgdGhlIHBsb3RzIG9mZgogIGdncGxvdChhZXMoeD1sb2cobGFtYmRhKSwgeT0gbWVhbmN2bSwgIHltaW49bG93LCB5bWF4PWhpZ2gpKSArCiAgZ2VvbV9yaWJib24oYWxwaGE9LjI1KSArCiAgZ2VvbV9saW5lKGFlcyhjb2xvcj1hcy5jaGFyYWN0ZXIoYWxwaGEpKSkgKwogIGZhY2V0X3dyYXAofiBhcy5jaGFyYWN0ZXIoYWxwaGEpKSArCiAgIGNvb3JkX2NhcnRlc2lhbih4bGltPShjKC01LDApKSkgKwogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQ9bG9nKGxhbWJkYS5taW4ubWVhbikpLCBhbHBoYT0uNSwgZGF0YT1tZXRfYW10X3N1bW1hcnlfbGFtYmRhKSArCiAgZ2VvbV92bGluZShhZXMoeGludGVyY2VwdD1sb2cobGFtYmRhLm1pbi5oaWdoKSksIGFscGhhPS41LCBkYXRhPW1ldF9hbXRfc3VtbWFyeV9sYW1iZGEsIGNvbG9yPSJibHVlIikgCgpgYGAKCgpNYWtlIGEgcGxvdCBvZiBNU0UgYXQgbWluaW11bSBsYW1iZGEgZm9yIGVhY2ggYWxwaGEKCmBgYHtyfQptZXRfYW10X3N1bW1hcnlfY3ZtICU+JSAKICBncm91cF9ieShhbHBoYSkgJT4lCiAgZmlsdGVyKHJhbmsobWVhbmN2bSwgdGllcy5tZXRob2QgPSAiZmlyc3QiKT09MSkgJT4lCiAgZ2dwbG90KGFlcyh4PWFscGhhLHk9bWVhbmN2bSx5bWluPWxvdyx5bWF4PWhpZ2gpKSArCiAgZ2VvbV9yaWJib24oY29sb3I9TkEsIGZpbGw9ImdyYXk4MCIpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludCgpCmBgYApub3QgYSBwYXJ0aWN1bGFyIGxhcmdlIGRpZmZlcmVuY2UgaGVyZSBhZnRlciAwLjIKClBsb3QgdGhlIG51bWJlciBvZiBuemVybyBjb2VmZmljaWVudHMKCmBgYHtyfQptZXRfYW10X3N1bW1hcnlfbGFtYmRhICU+JQogIHVubmVzdChjKGxhbWJkYSwgbnplcm8pKSAlPiUKICBncm91cF9ieShhbHBoYSkgJT4lCiAgZmlsdGVyKGFicyhsYW1iZGEubWluLm1lYW4tbGFtYmRhKT09bWluKGFicyhsYW1iZGEubWluLm1lYW4tbGFtYmRhKSkgICkgJT4lCiAgdW5ncm91cCgpICU+JQoKZ2dwbG90KGFlcyh4PWFzLmNoYXJhY3RlcihhbHBoYSksIHk9bnplcm8pKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZ3RpdGxlKCJOdW1iZXIgb2Ygbm9uLXplcm8gY29lZmZpY2VudHMgYXQgbWluaW11bSBsYW1iZGEiKSArCiAgeWxpbSgwLDM2KQpgYGAKT0sgbGV0J3MgZG8gcmVwZWF0ZWQgdGVzdCB0cmFpbiBzdGFydGluZyBmcm9tIHRoZXNlIENWIGxhbWJkYXMKCmBgYHtyfQptdWx0aV90dCA8LSBmdW5jdGlvbihsYW1iZGEsIGFscGhhLCBuPTEwMDAwLCBzYW1wbGVfc2l6ZT0zNiwgdHJhaW5fc2l6ZT0zMCwgeCwgeT1sZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCkgewogIHByaW50KGxhbWJkYSkKICBwcmludChhbHBoYSkKdHQgPC0KICB0aWJibGUocnVuPTE6bikgJT4lCiAgbXV0YXRlKHRyYWluPW1hcChydW4sIH4gc2FtcGxlKDE6c2FtcGxlX3NpemUsIHRyYWluX3NpemUpKSkgJT4lCiAgbXV0YXRlKGZpdD1tYXAodHJhaW4sIH4gZ2xtbmV0KHg9eFsuLF0sIHk9eVsuXSwgbGFtYmRhID0gbGFtYmRhLCBhbHBoYSA9IGFscGhhICkpKSAlPiUKICAKICBtdXRhdGUocHJlZD1tYXAyKGZpdCwgdHJhaW4sIH4gcHJlZGljdCgueCwgbmV3eCA9IHhbLS55LF0pKSkgJT4lCiAgbXV0YXRlKGNvcj1tYXAyX2RibChwcmVkLCB0cmFpbiwgfiBjb3IoLngsIHlbLS55XSkgICkpICU+JQogIG11dGF0ZShNU0U9bWFwMl9kYmwocHJlZCwgdHJhaW4sIH4gbWVhbigoeVstLnldIC0gLngpXjIpKSkgJT4lCiAgc3VtbWFyaXplKAogICAgbnVtX25hPXN1bShpcy5uYShjb3IpKSwgCiAgICBudW1fbHRfMD1zdW0oY29yPD0wLCBuYS5ybT1UUlVFKSwKICAgIGF2Z19jb3I9bWVhbihjb3IsIG5hLnJtPVRSVUUpLAogICAgYXZnX01TRT1tZWFuKE1TRSkpCnR0Cn0KCmFtdF9maXRfdGVzdF90cmFpbiA8LSBtZXRfYW10X3N1bW1hcnlfbGFtYmRhICU+JSAKICBzZWxlY3QoYWxwaGEsIGxhbWJkYS5taW4ubWVhbikKCmFtdF9maXRfdGVzdF90cmFpbiA8LSBtZXRfYW10X211bHRpQ1YgJT4lCiAgZmlsdGVyKHJ1bj09MSkgJT4lCiAgc2VsZWN0KGFscGhhLCBmaXQpICU+JQogIHJpZ2h0X2pvaW4oYW10X2ZpdF90ZXN0X3RyYWluKQoKYW10X2ZpdF90ZXN0X3RyYWluIDwtIGFtdF9maXRfdGVzdF90cmFpbiAlPiUKICBtdXRhdGUocHJlZF9mdWxsPW1hcDIoZml0LCBsYW1iZGEubWluLm1lYW4sIH4gcHJlZGljdCgueCwgcz0ueSwgbmV3eD1tZXRfYW10LlBDcykpLAogICAgICAgICBmdWxsX1I9bWFwX2RibChwcmVkX2Z1bGwsIH4gY29yKC54LCBsZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCkpLAogICAgICAgICBmdWxsX01TRT1tYXBfZGJsKHByZWRfZnVsbCwgfiBtZWFuKChsZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZC0ueCleMikpKSAlPiUKICAKICBtdXRhdGUodHQ9bWFwMihsYW1iZGEubWluLm1lYW4sIGFscGhhLCB+IG11bHRpX3R0KGxhbWJkYT0ueCwgYWxwaGE9LnksIHg9bWV0X2FtdC5QQ3MpKSkKCgoKKGFtdF9maXRfdGVzdF90cmFpbiA8LSBhbXRfZml0X3Rlc3RfdHJhaW4gJT4lIHVubmVzdCh0dCkpCmBgYAoKYGBge3J9CmFtdF9maXRfdGVzdF90cmFpbiAlPiUKICBnZ3Bsb3QoYWVzKHg9YWxwaGEpKSArCiAgZ2VvbV9saW5lKGFlcyh5PWF2Z19jb3IpLCBjb2xvcj0icmVkIikgKwogIGdlb21fcG9pbnQoYWVzKHk9YXZnX2NvciksIGNvbG9yPSJyZWQiKSArCiAgZ2VvbV9saW5lKGFlcyh5PWF2Z19NU0UpLCBjb2xvcj0iYmx1ZSIpICsKICBnZW9tX3BvaW50KGFlcyh5PWF2Z19NU0UpLCBjb2xvcj0iYmx1ZSIpCmBgYAphbHBoYSBvZiAwLjggdG8gMS4wIGFyZSB2ZXJ5IHNpbWlsYXIgYW5kIGFyZSB0aGUgYmVzdCBoZXJlLgoKIyMgbG9vayBhdCBmaXQ6CgpgYGB7cn0KYWxwaGFfYW10IDwtIC44CgpiZXN0X2FtdCA8LSBhbXRfZml0X3Rlc3RfdHJhaW4gJT4lIGZpbHRlcihhbHBoYSA9PSBhbHBoYV9hbXQpIApiZXN0X2FtdF9maXQgPC0gYmVzdF9hbXQkZml0W1sxXV0KYmVzdF9hbXRfbGFtYmRhIDwtIGJlc3RfYW10JGxhbWJkYS5taW4ubWVhbgoKYW10X2NvZWYudGIgPC0gY29lZihiZXN0X2FtdF9maXQsIHM9YmVzdF9hbXRfbGFtYmRhKSAlPiUgCiAgYXMubWF0cml4KCkgJT4lIGFzLmRhdGEuZnJhbWUoKSAlPiUgCiAgcm93bmFtZXNfdG9fY29sdW1uKHZhcj0iUEMiKSAlPiUKICByZW5hbWUoYmV0YT1gMWApCiAgCmFtdF9jb2VmLnRiICU+JSBmaWx0ZXIoYmV0YSE9MCkgJT4lIGFycmFuZ2UoYmV0YSkKCmBgYAoKcHJlZCBhbmQgb2JzCmBgYHtyfQpwbG90KGxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkLCBiZXN0X2FtdCRwcmVkX2Z1bGxbWzFdXSkKY29yLnRlc3QobGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQsIGJlc3RfYW10JHByZWRfZnVsbFtbMV1dKSAjLjczNgpiZXN0X2FtdCRmdWxsX01TRQpgYGAKCiMjIFBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkCgpgYGB7cn0KYW10X3ZhcnMgPC0gYW10X2NvZWYudGIgJT4lIAogIGZpbHRlcihiZXRhICE9MCwgUEMhPSIoSW50ZXJjZXB0KSIpICU+JQogIHB1bGwoUEMpICU+JSBjKCJsZWFmX2F2Z19zdGQiLCAuKQoKYW10X3JlbGltcCA8LSBsZWFmbGVuZ3RoICU+JSBzZWxlY3QobGVhZl9hdmdfc3RkKSAlPiUgY2JpbmQobWV0X2FtdC5QQ3MpICU+JSBhcy5kYXRhLmZyYW1lKCkgJT4lIGRwbHlyOjpzZWxlY3QoYWxsX29mKGFtdF92YXJzKSkgJT4lCiAgY2FsYy5yZWxpbXAoKSAKCmFtdF9jb2VmLnRiIDwtIGFtdF9yZWxpbXBAbG1nICU+JSBhcy5tYXRyaXgoKSAlPiUgYXMuZGF0YS5mcmFtZSgpICU+JQogIHJvd25hbWVzX3RvX2NvbHVtbigiUEMiKSAlPiUKICByZW5hbWUoUHJvcFZhcl9tZXRfYW10PVYxKSAlPiUKICBmdWxsX2pvaW4oYW10X2NvZWYudGIpICU+JQogIGFycmFuZ2UoZGVzYyhQcm9wVmFyX21ldF9hbXQpKQoKYW10X2NvZWYudGIKCmBgYAoKQ2hlY2tvdXQgdGhlIHJvdGF0aW9ucy4gIAoKYGBge3J9Cm1ldF9hbXRfcm90YXRpb25fb3V0IDwtIG1ldF9hbXQuUENfcm90YXRpb24gJT4lIAogIHBpdm90X2xvbmdlcigtbWV0YWJvbGl0ZSwgbmFtZXNfdG89IlBDIiwgdmFsdWVzX3RvPSJsb2FkaW5nIikgJT4lCiAgZmlsdGVyKFBDICVpbiUgZmlsdGVyKGFtdF9jb2VmLnRiLCBiZXRhIT0wKSRQQyApICU+JQogIGdyb3VwX2J5KFBDKSAlPiUKICAgIGZpbHRlcighc3RyX2RldGVjdChtZXRhYm9saXRlLCIuKihsZWFmfHJvb3QpX1swLTldKiQiKSkgJT4lCiAgZmlsdGVyKGFicyhsb2FkaW5nKSA+PSAwLjA1KSAlPiUKICBsZWZ0X2pvaW4oYW10X2NvZWYudGIsIGJ5PSJQQyIpICU+JQogIGFycmFuZ2UoZGVzYyhhYnMoYmV0YSkpLCBkZXNjKGFicyhsb2FkaW5nKSkpICU+JQogIG11dGF0ZShvcmdhbj1pZmVsc2Uoc3RyX2RldGVjdChtZXRhYm9saXRlLCAiX2xlYWZfIiksICJsZWFmIiwgInJvb3QiKSwKICAgICAgICAgdHJhbnNmb3JtYXRpb249InJhdyIsCiAgICAgICAgIG1ldGFib2xpdGU9c3RyX3JlbW92ZShtZXRhYm9saXRlLCAibWV0X2FtdF8ocm9vdHxsZWFmKV8iKSwKICAgICAgICAgbWV0YWJvbGl0ZV9lZmZlY3Rfb25fbGVhZj1pZmVsc2UoYmV0YSpsb2FkaW5nPjAsICJpbmNyZWFzZSIsICJkZWNyZWFzZSIpKQptZXRfYW10X3JvdGF0aW9uX291dCAlPiUgIHdyaXRlX2NzdigiLi4vb3V0cHV0L0xlYWZfYXNzb2NpYXRlZF9tZXRhYm9saXRlc19yYXcuY3N2IikKCm1ldF9hbXRfcm90YXRpb25fb3V0CmBgYAoK